Nonlinear multi-magnon scattering in artificial spin ice

Magnons, the quantum-mechanical fundamental excitations of magnetic solids, are bosons whose number does not need to be conserved in scattering processes. Microwave-induced parametric magnon processes, often called Suhl instabilities, have been believed to occur in magnetic thin films only, where quasi-continuous magnon bands exist. Here, we reveal the existence of such nonlinear magnon-magnon scattering processes and their coherence in ensembles of magnetic nanostructures known as artificial spin ice. We find that these systems exhibit effective scattering processes akin to those observed in continuous magnetic thin films. We utilize a combined microwave and microfocused Brillouin light scattering measurement approach to investigate the evolution of their modes. Scattering events occur between resonance frequencies that are determined by each nanomagnet’s mode volume and profile. Comparison with numerical simulations reveals that frequency doubling is enabled by exciting a subset of nanomagnets that, in turn, act as nanosized antennas, an effect that is akin to scattering in continuous films. Moreover, our results suggest that tunable directional scattering is possible in these structures.


SUPPLEMENTARY NOTE 1. SCANNING ELECTRON MICROSCOPY
Supplementary Figure 1 shows scanning electron microscopy images of the ASI under study. The size discrepancy between the design and the fabricated nanoelements, as well as variability between elements, is smaller than 2 nm.

SUPPLEMENTARY NOTE 3. POWER THRESHOLD FITS
We approximate the threshold behavior using the auto-oscillator theory [1,2] to estimate the threshold power. The theory predicts that a generalized auto-oscillator has a power described by where ξ = P ω /P th is called the supercriticallity parameter and it is a ratio of the microwave power to the microwave threshold power. Because this equation predicts an output power that increases asymptotically, it is difficult to fit a threshold, especially for experimental data with noise. A better representation consists in plotting the microwave power, so that the resulting data can be fitted with a line. We found that for this data set, smaller error bars in the fitted quantities were obtained by setting zero as the noise-floor, < P >.
A particular difficulty with our data set is the lack of points after the threshold is detected.
For this reason, we consider an increasing number of data points and determine the threshold based on the best fit. An example of the best fit is shown in Supplementary Fig. 2. The threshold is then determined as the point where the fit crosses < P > −P = 0. We obtain a threshold power of 100 mW ± 3 mW (20 dBm ± 5 dBm) from this fit.

SUPPLEMENTARY NOTE 4. CHIRPED SIMULATIONS
The main idea behind chirped micromagnetic simulations is to save computational time.
Experimentally, the spectra are obtained by capturing a time-averaged spectrum for each experimental condition: bias field, microwave frequency, and microwave power. The signalto-noise ratio is improved by averaging multiple spectra. This technique can be reproduced numerically but at a high computational cost. For example, suppose a microwave excitation at 1 GHz and a requirement to obtain a spectrum with a frequency resolution of 100 MHz, comparable to BLS. According to the Nyquist criterion, the minimum sampling time will be 0.5 ns from the microwave frequency. However, the full spectrum up to 20 GHz, e.g., in Fig. 2, requires a sampling time of at least 25 ps, increasing the trace to 400 points.
Considering that a single oscillation of the microwave frequency is insufficient to achieve a steady-state, a typical stabilization time of 10 ns must be considered. This already introduces an 11-ns-long simulation for each experimental condition. Finally, we note that the sampling time is different from the time step of the partial differential equation method solver used. This time step is typically an order of magnitude smaller, i.e., 1 ps. Overall, at least 11,000 iterations would be required per experimental condition. This argumentation illustrates the enormous computational resources required to reproduce Fig. 2 by simulating the experimental method.
An alternative solution is to use a smooth change in the external conditions to approximately recover the experimental results. This technique economizes the use of resources with acceptable accuracy, as clearly illustrated in Ref. [3]. Here, we introduce this technique for the case of microwave excitation in artificial spin ices, where the rate of the chirp was adjusted to minimize transient effects. We also perform the simulation in blocks due to the large amount of data required to compute the smoothed pseudo-Wigner-Ville distribution [4,5]. As such, we can perform only two 400 ns long simulations to explore a range of 20 GHz of microwave excitations. Brute-force simulations replicating the BLS experimental approach would instead require a net simulated time of 18 µs, almost 46 times longer simulation time.
As mentioned above, this technique is prone to artifacts. An example is the broad spectrum in Fig. 5d at about 6 GHz. This feature appears because of a transition in the oscillations and thus the transient dynamics towards stabilization are captured as well. Despite these artifacts, most features are accurately described, allowing us to choose specific parameters to run dedicated simulations and resolve the magnetic mode volumes.

SUPPLEMENTARY NOTE 5. ASI SUBLATTICE MODES
We have previously shown that each branch in the thermal modes can be related to each sublattice, e.g., [6][7][8]. This is confirmed in the ASI sublattice mode simulations shown in Supplementary Fig. 3. We have argued in the manuscript that parametric pumping/ higher harmonic generation are achieved in our experiments due to the strong microwave excitation. In addition, we have found that a frequency-doubled mode is excited in vertical nanoislands that are not due to higher harmonic generation, nor such modes are observed in ASI sublattices. The frequencydoubled mode is argued to appear as a consequence of higher harmonic generation in the horizontal nanoislands and these, in turn, acting as nanoantennas that excite the vertical nanoislands at the double frequency. Here, we show that indeed the horizontal nanoislands excite the vertical elements due to higher harmonic generation.
A tell-tale of higher harmonic generation is ellipticity in the magnetization precession array. This demonstrates that the magnetic state is primarily determined by the internal energy, in particular shape anisotropy. The thermal spectrum as a function of field shown in Supplementary Fig. 5b, respectively, is also similar to those obtained for the horizontal and vertical sublattices. So far, these results indicate that the salient features of the magnon spectrum are only slightly modified by the magnetostatic coupling between nanoislands.
This is well known in the literature, and only edge modes are expected to interact between nanoislands [9].
The modes obtained under microwave excitation are shown in Fig. 5c for f = 4 GHz at B ω = 5 mT and Supplementary Fig. 5d for f = 19 GHz at B ω = 50 mT. Comparison with the results of ASI sublattices in the main text, Fig. 7, we observe essentially the same modes, barring minor mode profile distortion in the ASI sublattice due to the magnetostatic field.
With these simulations, we further corroborate that parametric pumping occurs within the nanoislands and that the mode profiles are mainly determined by the internal field in the nanoislands.

SUPPLEMENTARY NOTE 8. 3D SIMULATIONS
To investigate the impact of our cell dimensions, we performed more detailed simulations using a reduced cell size of dimensions 1 nm ×1 nm ×1 nm. These cells have three times smaller lateral dimensions compared to those used in Fig. 5 and are also cubic, which improves the accuracy of the computation. The total number of cells for our geometry is 690 × 690 × 15 resulting in ≈ 7.1 million cells. In contrast, the simulations presented in the main text require only ≈ 53 thousand cells. The significant increase in total cells results in a much longer simulation time. For example, a single microwave excitation to obtain mode profiles required 72 hours compared to 2 hours for the results presented in the main text.
In Supplementary Fig. 6, we present results from simulations with the reduced cell size for both the thermal modes and the case of a microwave excitation of f = 4 GHz and B ω = 5 mT. The mode profiles are shown in panels a and b for the m y component. We observe a similar spatial distribution as in the main text Fig. 5b and c. However, we find a quantitative discrepancy in the modes' frequency. In particular, the edge mode is found at 8.08 GHz and the bulk mode at 15 GHz, cf. 7.89 GHz and 13.75 GHz, respectively, shown in Figure 5 in the main text. The discrepancy originates from a slightly different distribution of the magnetization state at that particular field. The nanomagnets find themselves in a preferentially horizontal configuration rather than a tilted configuration for the 3D simulations presented here. This is due to the rougher edges and lack of volume resolution in the results with larger cells.  we find that the magnetization in the horizontal nanoislands is not excited in this simulation.
Instead, the vertical nanoislands can more easily couple to the microwave field due to the reduced anisotropy in the y-direction. We find, nonetheless, a qualitative agreement in the excitation of the direct 4 GHz mode and the appearance of a 2f harmonic that is spatially located in a different section of the sample.